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STEADY-STATE SIMULATION OF REFLECTED 
BROWNIAN MOTION AND RELATED STOCHASTIC 

NETWORKS 

By Jose Blanchet and Xinyun Chen 

Columbia University 

This paper develops the first class of algorithms that enables un- 
biased estimation of steady-state expectations for multidimensional 
reflected Brownian motion. In order to explain our ideas we first con- 
sider the case of compound Poisson input. In this case, we analyze 
the complexity of our procedure as the dimension of the network in- 
creases and show that under certain assumptions the algorithm has 
polynomial expected termination time as the dimension increases. 
Our methodology includes procedures that are of interest beyond 
steady-state simulation and reflected processes. For instance, we use 
wavelets to construct a piece-wise linear function that can be guar- 
anteed to be within e distance (deterministic) in the uniform norm 
to Brownian motion in any compact time interval. 

1. Introduction. This paper studies simulation methodology that al- 
lows to estimate, without any bias, steady-state expectations of reflected 
(also known as constrained) stochastic networks. Our algorithms are pre- 
sented with companion rates of convergence for the steady-state computa- 
tions of reflected stochastic networks. Reflected stochastic networks, as we 
shall explain, are very important for the analysis of queueing systems. How- 
ever, in order to motivate the models that we study let us quickly review a 
formulation introduced by Kella [1996]. 

Consider a network of d queueing stations indexed by {1, 2, d}. Suppose 
that jobs arrive to the network according to a Poisson process with rate A, 
denoted by {N (t) : t > 0). Specifically, the k-th arrival brings a vector of job 
requirements W (k) = {Wi (fc) , Wd {k))^ which are non-negative random 
variables (r.v.'s) and that add to the workload at each station right at the 
moment of arrival. So, if the fc-th arrival occurs at time t, the workload of 
the i-th station (for i € {1, d}) increases by Wi {k) units right at time t. 
We assume that W = (W (/c) : A; > 1) is a sequence of i.i.d. (independent 
and identically distributed) non-negative r.v.'s. For fixed fc, the components 
of W (fc) are not necessarily independent, however, W is assumed to be 
independent of (•). 
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Throughout the paper we shall use boldface notations to write vector 
quantities, which are encoded as columns. For instance, we write y = (yi, ■■.,yd) 

The total amount of work that arrives to the i-th. station by time t is 
denoted by 

N{t) 

j.it) = Y.w,ik). 

k=l 

Now, assume that the workload at the i-th station is processed as a fluid by 
the server at a rate rj, continuously in time. This means that if the workload 
in the i-th station remains strictly positive during the time interval \t,t + dt] 
then the output from station i during this time interval equals Vidt. In 
addition, suppose that a proportion j > of the fluid processed by the i-th 
station is circulated to the j-th server. We have that Qij < 1, Qi^i = 0, 

and we define Qi^ = 1 — Yl'j=i QiJ- "^^^ proportion (5i,o corresponds to the 
fluid that goes out of the network from station i. 

The dynamics stated in the previous paragraph are expressed formally in 
differential notation as follows. Let Yi (t) denote the workload content of the 
i-th station at time t, then given Yi (0), we write 

dY, (t) = dJi (t) - nl {Y, (t) >Q)dt+Y^ Qj^i^jl (^j (0 > 0) dt 
= dJi {t) - ridt + ^ Qj,irjdt 

(1) + m {Yi (t) = 0)dt-Y, Qi,^^il (^i (*) = 0) dt 

for i G It is immediate that the resulting vector- valued workload 

process, Y {t) = (Yi {t) , ...,1^ (i))^, is Markovian. The differential equation 
as defined by (1) admits a unique right-continuous-with-left-limits (ROLL) 
piecewise linear solution. This can be easily established by elementary meth- 
ods and we shall comment on far reaching extensions shortly. 

The equations given in (1) take a neat form in matrix notation. This 
notation is convenient to discuss stability issues and other topics which are 
related to the steady-state simulation problem that concerns this paper. In 
particular, if we let r = (ri, r^) be the column vector corresponding to 
the service rates, write R = {I — Q)'^ , and define 

X (t) = J (t) - Rrt, 

where J (t) is a column vector with its i-th component equal to Jj (t), then 
an equivalent form of equation (1) is expressed as 

(2) Y(t) = Y(0) + X(t) + i2L(t), 
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where L (t) is a column vector with its i-th component equal to 

Li (t) = [ rj {Yi {s) = 0) ds. 
Jo 

As mentioned earlier Y = (Y (t) : t > 0) is a Markov process, to ensure 
stability let us assume that — t- as n — t- oo. This assumption in particular 
implies that the network is open in the sense that for each i such that Aj > 
there exists a path {ii,i2, ■.■,ik) with ik = 0, k < d and satisfying that 
KQi,iiQii,i2---Qik-i,ik > 0. In addition, under this assumption the matrix 
exists and it has non- negative components. Moreover, suppose that 
i?~^ii^X(l) < - inequalities involving vectors are understood component- 
by-component throughout the paper. It follows from results of Kella [1996] 
that Y (t) =^ Y (oo) as t — 7- oo, where Y (oo) is a r.v. with the (unique) 
stationary distribution of Y (•). 

Our first contribution in this paper is to develop an exact sampling algo- 
rithm (i.e. simulation without bias) for Y (oo). This algorithm is developed 
in Section 2 of this paper under the assumption that W (k) has a finite mo- 
ment generating function. In addition, we analyze the rate of convergence 
(measured in terms of expected random numbers generated) to terminate the 
algorithm as d increases. 

The workload process (Y (t) : t > 0) is a particular case of a reflected 
(or constrained) stochastic network. Although the models introduced in the 
previous paragraphs are interesting in their own right, our main interest is 
on the study of steady-state simulation techniques for reflected Brownian 
motion, which is obtained by abstracting the construction formulated in 
(2). This abstraction is presented in terms of a Skorokhod problem, which 
we now describe. Given a matrix R, and a general ROLL input process 
X = (X (t) :t>0), which we eventually shall take as a Brownian motion 
with given drift v = i?X(l) and covariance matrix S = Var(X(l)), the 
Skorokhod problem consists in finding a pair of processes (Y,L) satisfying 
precisely equation (2), subject to the initial condition Y (0) > and the 
following, 



Skorokhod Problem Constraints: 

i) Y (t) > for each t, 
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ii) Li (•) non-decreasing for each i G {1, ...d} and Lj (0) = 0, 

iii) Jq Yi (s) dLi (s) = for each t. 

It turns out that if the input process X(-) is Brownian motion and if 
the matrix i? is a so-cahed M-matrix, which is the same as saying that 
exists and it has non-negative components, then there exists a strong 
solution (i.e. path-by-path and not only in law) to the stochastic differential 
equation (SDE) (2) subject to the Skorokhod Problem Constraints i) to iii), 
and the initial condition Y (0) . This was proved by Harrison and Rieman 
[1981], who introduced the notion of reflected Brownian motion (RBM). In 
this paper we assume that R is an M-matrix. In this setting, it follows, 
see for instance, Kella and Whitt [1996] that R~^^ < is a necessary and 
sufficient condition for stability; see also Chapter 7 of Chen and Yao [2001] 
and references therein. 

Dupuis and Ishii [1993] extended the class of input processes for which a 
strong solution to the Skorokhod problem can be shown to exist. Existence 
and uniqueness of the solution to the Skorokhod problem is an area of current 
research motivated by applications both in Operations Research and other 
areas, see for instance Dupuis and Ramanan [2000], Bhardwaj and Williams 
[2009]. 

In the setting of RBM, Harrison and Rieman [1981] were motivated by 
the fact that in great generality (only under the existence of variances of 
service times and inter-arrival times) so-called generalized Jackson networks 
(which are single-server queues connected with Markovian routing) converge 
weakly to reflected Brownian motion in a heavy traffic asymptotic envi- 
ronment, see Reiman [1984]. Therefore, reflected Brownian motion (RBM) 
plays a central role in queueing theory. Moreover, recent papers, for instance 
Gamarnik and Zeevi [2006] and Budhiraja and Lee [2009], have shown that 
convergence occurs also at the level of steady-state distributions. 

Our second and most important contribution in this paper is the devel- 
opment of an algorithm that allows to estimate with no bias Eg (Y (oo)) 
for positive and continuous functions g (•). Moreover, given e > 0, we pro- 
vide a simulation algorithm that outputs a random variable Y^ (oo) that 
can be guaranteed to be within e distance (say in the Euclidian norm) from 
an unbiased sample Y (oo) from the steady-state distribution of RBM. This 
contribution is developed in Section 3 of this paper. We show that the num- 
ber of Gaussian random variables generated to produce 1^ (oo) is of order 
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0(e~'^'^~^ log(l/e)) as e \ 0, where ac is a constant only depending on the 
covariance matrix of the Brownian motion (see Section 3.4)- In the special 
case when the d-Brownian motion has non-negative correlations, the number 
of random variables generated is of order 0{e~'^~^ \og{l/£)). 

Our methods allow to estimate without bias Eg (Y [ti) , Y (t2) , Y (tm.)) 
for a positive function g (•) continuous almost everywhere and for any < 
ti < t2 < ... < tra- Simulation of RBM has been studied in the literature. 
In the one dimensional setting it is not difficult to sample RBM exactly; 
this follows, for instance, from the methods in Devroye [2009]. The paper of 
Asmussen et al. [1995] also studies the one dimensional case and provides an 
enhanced Euler-type scheme with an improved convergence rate. The work 
of Burdzy and Chen [2006] provides approximations of reflected Brownian 
motion with orthogonal reflection (the case in which R = I). 

On the side of steady-state estimation, the work of Dai and Harrison 
[1992] provides numerical methods for approximating the steady-state ex- 
pectation by numerically evaluating the density of Y(oo). In contrast to 
our methods, Dai and Harrison's procedure is based on finite elements and 
thus it is non-randomized, which is in some sense preferable. However, the 
theoretical justification of Dai and Harrison's algorithm relies on a conjec- 
ture that is believed to be true but has not been rigorously established, see 
Dai and Dicker [2011]. In addition, no rate of convergence is given for this 
procedure, even assuming that the conjecture is true. 

Finally, we briefly discuss some features of our procedure and our strategy 
at a high level. There are two sources of bias that arise in the setting of 
steady-state simulation of RBM. First, simulating the process Y requires 
implicitly the simulation of the local-time-like process L. The fact that the 
reflection matrix R is not the identity makes this task very difficult. This 
issue is present even in finite time horizon. The second issue is, naturally, that 
we are concerned with steady-state expectations which inherently involve, 
in principle, an infinite time horizon. 

In the reflected compound Poisson case we can simulate the solution of 
the Skorokhod problem in any finite interval exactly and without any bias. 
So we can concentrate on the problem of removing the bias issues arising 
from the infinite horizon, i.e. we can concentrate on the initial transient 
problem. This is actually our main motivation for considering these types of 
models. 
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The strategy that we pursue is based on Dominated Couphng From The 
Past (DCFTP). This technique was proposed by Kendall [2004], following 
the introduction of Coupling From The Past by Propp and Wilson [1996]. 
The idea behind DCFTP is to construct a suitable dominating upper and 
lower bound processes that can be simulated in stationarity and backwards 
in time. We take the lower bound to be the process identically equal to zero. 
We use results from Kella and Whitt [1996], see also Ramasubramanian 
[2000], to construct an upper bound process based on the solution of the 
Skorokhod problem with reflection matrix R = I. It turns out that simu- 
lation of the stationary upper-bound process backwards involves sampling 
the infinite horizon maximum (component-by-component) from t to infinity 
of a (i-dimensional compound Process with negative drift. We use sequential 
acceptance / rejection techniques (based on a exponential tilting distribu- 
tions used in rare-event simulation) to simulate from such infinite horizon 
maximum process. 

Then we turn to RBM. A problem that arises, in addition to the discretiza- 
tion error given the continuous nature of Brownian motion, is the fact that 
in dimensions higher than one (as in our setting) RBM never reaches the 
origin. Nevertheless, it will be arbitrarily close to the origin and we shall 
certainly leverage off this property to obtain simulation that is guaranteed 
to be e-close to a genuine steady-state sample. Now, in order to deal with 
the discretization error we use wavelet-based techniques. We take advantage 
of a well-known wavelet construction of Brownian motion, see Steele [2001]. 

Instead of simply simulating Brownian motion using the wavelets, which 
is the standard practice, we simulate the wavelet coefficients jointly with 
suitably defined stopping times. Consequently, we are able to guarantee that 
our wavelet approximation is e-close (note that e is deterministic and defined 
by the user) in the uniform metric to Brownian motion in any compact time 
interval (see Section 3.2). 

Finally, we use the fact that the process Y arising as the solution to the 
Skorokhod problem, as a function of the input process X, is Lipschitz con- 
tinuous with a computable Lipschitz constant, under the uniform topology. 
These observations together, combined with an additional randomization, in 
the spirit of Beskos et al. [2012] allow to estimate expectation with no bias. 

We strongly believe that the use of tolerance-enforced coupling based on 
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wavelet constructions, as we illustrate here, can be extended more broadly 
in the numerical analysis of the Skorokhod and related problems. 

Numerical experiments are pursued in a companion paper in which we also 
discuss further implementation issues of importance. We prefer to confine 
the scope of this paper only to the conceptual and theoretical properties 
behind our algorithm. 

The rest of the paper is organized as follows. In Section 2, we consider 
the problem of exact simulation from the steady-state distribution of the 
reflected compound Process discussed earlier. This section allows to explain 
the main strategy to be used for the reflected Brownian motion case, which 
is studied in Section 3. 

2. Exact Simulation of Reflected Compound Poisson Processes. 

The model that we consider has been explained at the beginning of the 
Introduction. We summarize the assumptions that we shall impose next. 

Assumptions: 

Al) The matrix R is invertible and has non-negative components, 
A2) R~^E'K (1) < (recall that inequalities apply component-by-component 
for vectors), 

A3) There exists 6 >0, 6 eR"^ such that 

logS[exp (^^W(fc))] < 0. 



We have commented on Al) and A2) in the Introduction. Assumption 
A3) is important in order to simulate from the steady-state upper bound 
process that we shall construct to apply DCFTP. 

In addition to Al) to A3) we shall assume that one can simulate from 
exponential tilting distributions associated to the marginals of W(/c). That 
is, we can simulate from Pg. (•) such that 

Pe, {Wi{k)edyi,...,Wd{k)€dy^) 
^exp [diWi [k)) 

where 0j G M and E exp {9iWi (k)) < oo. This is not a strong assumption 
provided one can simulate from W (k) directly. 
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Let us briefly explain our program, which is based on DCFTP. First, we 
will construct a stationary dominating process (Y+ (s) : — oo < s < 0). The 
dominating process will be constructed coupled with a stationary process 
(Y (s) : — oo < s < 0), which satisfies the Skorokhod problem (2) and also 

(3) R^^Y {s) < R-^Y+ {s) , 

for each s < 0. We simulate the process Y"*" (•) up to a time — r < such that 
Y+ (— r) = 0. Following the tradition of the CFTP literature, we call a time 
— T such that Y^ (—t) = a coalescence time. Since Y (s) > 0, inequality 

(3) yields that Y (— r) = 0. The next and final step in our strategy is to 
evolve the solution of the Skorokhod problem (2) for r units of time forward, 
starting with zero initial condition, and using the same input that drives the 
construction of (Y (s) : — r < s < 0). The output is therefore Y (0) which 
is, of course, stationary. The precise algorithm will be summarized in Section 
2.2. 

So, a crucial part of the whole plan is the construction of Y"^ (•) together 
with a coupling that guarantees inequality (3). In addition, the coupling 
must be such that one can use the driving randomness that defines Y"*" (•) 
directly as an input to the Skorokhod problem (2) that is then used to evolve 
Y"*" (•). We shall first start by constructing a time reversed stationary version 
of a suitable dominating process Y"*". 

2.1. Construction of the Dominating Process. In order to construct the 
dominating process Y"*" (•) we first need the following result due to Kella and Whitt 
[1996]. 

Lemma 1 (Kella and Whitt '96). There exists z such that i?X(l) < z, 
R~^z < and such that if 

Z{t) =X(t) -zf, 
and ifY^ (•) satisfies the Skorokhod problem 

(4) dY+ (t) = dZ (t) + dL+ (t) , Y+ (0) = yo, 

Y+ (t) > 0, Y+ (t) dL+ (t) > 0, L+ (0) = 0, dLj (t) > 0. 

Then, < ii^^Y (t) < R~^Y+ (i) for all t > where Y (•) solves the 
Skorokhod problem 

dY (t) = dX (t) + RdL (t) , Y (0) = yo 
Y (t) > 0, Yj (t) dLj (t) > 0, Lj (0) = 0, dLj (t) > 0. 
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We note that computing z from the previous lemma is not difficult, one 
can simply pick z = E'X (1) + (51, where 1 = (1, 1)"^ and with 6 chosen 
so that < 6R~^1 < — i?~^£'X (1). In what follows we shall assume that z 
has been selected in this form. Note now that with such a selection we have 
that EZ (1) < 0. 

The Skorokhod problem corresponding to the dominating process can 
be solved explicitly. Indeed, it is not difficult to verify, see for instance 
Harrison and Rieman [1981], that if Y"^ (0) = 0, the solution of the Sko- 
rokhod problem (4) is given by 

Y+ (t) = Z (t) - min Z (n) = max (Z (t) - Z (n)), 

0<u<t 0<u<t 

where the running maximum is obtained component-by-component. 

In order to construct a stationary version of (•) backwards in time 
let us define a time-reversed version of Z(-). First, we extend Z(-) to a 
two-sided compound Poisson process, with Z (0) = 0. Then, because Z (•) 
has stationary and independent increments it follows that Z^ (•) defined via 
Z'*~ (t) = — Z (— t) has also stationary and independent increments. In fact, 
Z^ (•) and Z (•) are identical in distribution. 

Now, given T < define Z^ via Z|r (t) = {T + 1) for t > 0. In 
addition, given y > let Y^ (•, y) be the solution to the Skorokhod problem 

dY+ (t, y) = dZ^ (t) + dL+ {t, y) , Y+ (0, y) = y, 
Y+ {t,y) > 0, Y+ . (t,y) dL+ . (t,y) > 0, L+ . (0,y) = 0, dL+. (t,y) > 0. 

As noted earlier, we have that 

(5) Y+ {t, 0) = max (Zjr (t) - Z^ (n)). 

U<u<i 

The process Y"*" satisfying the Skorokhod problem with orthogonal reflection 
(it! = /) possesses a unique stationary distribution. So, we can construct a 
stationary version of (Y"*" (s) : —oo < s < 0) via 

(6) Y+(s)= lim Y+(-r + s,0). 

T-i.-oo 

The next proposition provides an explicitly evaluation of the limit in (6). 
Proposition 1. Given any t>0, 

(7) Y+(-t) = -Z{t) + max Z(n), 

t<u<oo 
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Proof. Expression (5) together with the definition of (•) yields, for 
T, s<0: 

Y+ (-r + s, 0) = max (Z^ (T + t) - Z^ (T + u) 

0<u<t 

= max (Z^(s)-Z^(T + ^.)) 

U<n<— T+s 

= max (Z^(s)-Z^(r + n)). 

T<u+T<s 

Now we let r = n + T, substitute Z"*~ (s) = — Z (— s) and obtain 
Y+ (-T + s, 0) = max (Z^ (s) - Z^ (r)) 

^ T<r<s 

= max(-Z(-s) + Z(-r)) 

T<r<s 

= -Z (-s) + max Z (-r) . 
r<r<s 

Letting — s = t > and — r = u > we obtain 

Y^(-T-i,0) = -Z(t) + max Z (n) . 

i<«<— T 

Now send —T — )• oo and arrive at (7), thereby obtaining the result. □ 



2.2. The Structure of the Main Simulation Procedure. We now are ready 
to explain our main algorithm to simulate unbiased samples from the steady- 
state distribution of Y. For this purpose, let us first define 

M{t) = max Z{u), 

t<u<oo 

for t > 0. Since E[Z (1)] < it follows that (M {t) : t > 0) is a well defined 
stochastic process. Let us for the moment assume that we can simulate M (•) 
on any compact time interval. We shall explain how to achieve this after we 
summarize our strategy in the next algorithm. 



Algorithm 1: Exact Sampling of Y (oo) 

Step 1: Let r > be any time for which M (r) = Z (r) and simulate 
(Z : < t < r). 
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Step 2: Set Xt:^ (t) = -Z (r - t) - z x (r - t) and compute (r, 0) 
which is obtained by evolving the solution Y_t- (•, 0) to the Skorokhod prob- 
lem 

dV-r (t, 0) = dXt; (t) + Rdl,_r (t, 0) 

(t, 0) > 0, Y^r,j {t, 0) dL_r,j {t, 0) > 0, L_^j (0, 0) = 0, dL_r,j {t, 0) > 0, 

for T units of time. 

Step 3: Output Y_,- (r, 0) which has the distribution of Y (oo). 

The time — r precisely corresponds to a coalescence time. The most inter- 
esting part that remains to discuss is the simulation of the stochastic object 
(M (t) : < t < t). Step 2 is straightforward to implement because the pro- 
cess Xt"^ (■) is piecewise linear and therefore the solution to the Skorokhod 
problem, namely Y_t- (-,0) is also piecewise linear. The gradients are sim- 
ply obtained by solving a sequence of linear system of equations which are 
dictated by evolving the ordinary differential equations given in (1). The 
following proposition summarizes the validity of this algorithm. 

Proposition 2. The previous algorithm terminates with probability one 
and its output is an unbiased sample from the distribution ofY (oo). 

Proof. It follows from the results in Kella [1996] that the distribution 
of Y"^ (oo) has an atom at the origin and therefore r < oo almost surely. 
Let T < and note that thanks to Lemma 1, for t S (0, |T|] 

R-^YT{t,0) < R-^Y+{t,0). 

In addition, by monotonicity of the solution to the Skorokhod problem in 
terms of its initial condition, see Ramasubramanian [2000], we also have 
(using the definition of Y+(T) from (6)) that 

Y+(t,o) < Y+(i,Y+(r)) = Y+(r + i). 

Therefore, 

Y+(r + t) = implies Yxit, 0) = 0. 
Consequently, if — T > r > 

Yt(|T|-t,0) = 0, 
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which in particular yields that Yq"(— T, 0) = Y_T-(r, 0). We then obtain that 
^hm YT(-r,0) = Y_,(r,0), 

1 — ^ — oo 

thereby concluding that Yt-(— r, 0) follows the distribution Y (oo) as claimed. 

□ 

2.3. Simulation of the Stationary Dominating Process. In order to sim- 
plify the explanation of the simulation procedure to sample (M (t) : t > 0) 
we introduce the following assumption. 

Assumption: A3b) Suppose that for every i there exists 9* £ (0,oo) 
such that 

{9*)^logEoexp{6*Z,{t)) = l. 

As customary, we use the notation Eq {■) or Pq {■) to indicate the conditioning 
Z(0) = 0. 

This assumption is a strengthening of Assumption A3) and it is known 
as Cramer's condition in the large deviations literature. As we shall explain 
later, it is possible to dispense this assumption and only work under As- 
sumption A3). For the moment, we continue under Assumption A3b). 

We wish to simulate (Z (t) : < t < r) where r is a time such that 
Z(r) = M(r) = maxZ(s). 

S>T 

Recall that — r is precisely the coalescence time. We also keep in mind that 
our formulation at the beginning of the Introduction implies that 

N{t) 

J (t) - i?rt - Zt = ^ W {k) - Rrt - zt. 

k=l 

Define 

fj, = Rr + z 

and let ;Uj < be the i-th component of fi. In addition, suppose that a 
constant m > is chosen to satisfy 

Condition 1: 

d 

(8) ^exp(-6'im) < 1. 

1=1 
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Now, define 

(9) = inf{t > : > m, for some i}. 

Using Condition 1 we can then propose the following procedure to simu- 
late r. The execution of the steps inside this procedure will be discussed 
momentarily. 

Algorithm 1.1: Simulating the Coalescence Time 

The output of this algorithm is (Z (i) : < t < r), and the coalescence 
time r. 

Set m according to (8). 

1. Set r = 0, Z (0) = 0, 

2. Generate an inter-arrival time U distributed Exp(A) and sample W = 
(VFi, Wd) independent of U. 

3. Let Z (r + = Z (r) - i/i for < t < [/ and Z (r + [/) = Z (r) + W - 
U (1. Then, reset r < — t + U . 

4. If there exists a component i, such that Wi — Ufn > —m, then return to 
Step 2. Otherwise, sample a Bernoulli I with parameter p = PQ(Tm < 
oo). 

5. If I = 1, simulate a new conditional path (C (t) : < t < T^) following 
the conditional distribution of (Z (t) : < t < Tm) given that Tm < oo 
and Z (0) = 0. Let Z (r + t) = Z (r) + C {t) for < f < T„ and reset 
T < — T+ Tm- Return to Step 2. 

6. Else, if / = 0, stop and return r along with the feed-in path (Z(t) : < t < r). 

We shall now explain how to execute the key steps in the previous algo- 
rithm, namely, Steps 4 and 5. 

2.3.1. Simulating a Path Conditional on Reaching a Positive Level in Fi- 
nite Time. The procedure that we shall explain now is an extension of the 
one dimensional procedure given in Blanchet and Sigman [2011]; see also the 
related one dimensional procedure by Ensor and Glynn [2000] . The strategy 
is to use acceptance / rejection. The proposal distribution is based on im- 
portance sampling by means of exponential tilting. In order to describe our 
strategy we need to introduce some notation. 



We think of the probability measure i-b (•) as defined on the canonical 
space of right-continuous with left-limits M"^- valued functions, namely, the 
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ambient space of (Z (t) : t > 0) which we denote hj Q = -D[o,oo) {^'^) ■ We en- 
dow the probabihty space with the Borel cj-field generated by the Skorokhod 
Ji topology, see Bilhngsley [1999]. We introduce our proposal distribution, 
Pq (•), defined on the space = -D[o,oo) (l^'^) x {1)2, ...,d}. We endow the 
probability space with the product cj-field induced by the Borel cr-field gen- 
erated by the Skorokhod Ji topology and all the subsets of {1, 2, d}. So, a 
typical element uj' sampled under Pq (•) is of the form a;'=((Z {t) : t > 0),K), 
where K G {1,2, ...,d}. The distribution of u' induced by Pq (•) is described 
as follows, first, 

(10) PHK = k)^^,:= /''P'^"" . . 

E-.iexp(-e'm) 

Now, given K = k, for every set A S (t(Z (s) : < s < t), 

P^ {A\K = k) = Eo (exp {eiZk (t)) I a) • 

So, in particular, the Radon-Nikodym derivative (i.e. the likelihood ratio) 
between the distribution of a; = (Z (s) : < s < t) under Pq (•) and Pq (•) is 
given by 

jpt <i 



The distribution of (Z (s) : s > 0) under Pq (•) is precisely the proposal 
distribution that we shall use to apply acceptance / rejection. It is straightfor- 
ward to simulate under Pq (•). First, sample K according to the distribution 

(10) . Then, conditional on K, the process Z (•) also follows a compound 
Poisson process. Indeed, given K = k, under Pq (•) it follows that J (t) can 
be represented as 

N'{t) 

(11) J(t)=J]W'(f), 

i=l 

where N' {■) is a Poisson process with rate A(/) (0, 0, 0^,, 0, 0), where 

(Pim, ■■■,Vd) = Pexp {riiWi + ... + 7]dWd) ■ 
In addition, we have that ii rj = {rji, ...,r/^)^, then 

(12) E (exp (, W (.)) \K = k)= ^(o^...^o^,.^o^...^o) • 



imsart-aap ver. 2011/11/15 file: Exact_Sampling_appl.tex date: February 10, 2012 



STEADY-STATE SIMULATION OF REFLECTED BROWNIAN MOTION 15 



So, conditional on K = k, we simply let 

N'{t) 

(13) Z{t)= J] W'(i)-//t. 

i=l 



Now, note that we can write 

d 



E', {Zk (t)) = EoiZk exp {OlZk {K = k) 

k=0 

^ {0*k) ^„ 

k=0 

where the last inequality follows a.s. by convexity of ip^ (■) and by definition 
of 0^. So, we have that (t) oo as t oo under Pq (•) by the Law of 
Large Numbers. Consequently Tm < oo a.s. under Pq (•). 

Finally, in order to show that we indeed can apply acceptance / rejection 
using the distribution of (Z (t) : t > 0) under P' (•), let us write Pq (•) to 
denote the conditional law of (Z (t) : < t < Tm) given that Tm < oo and 
Z (0) = 0. Then, we have that 

^(Z(t):0<t<r^) 
1 dP 

= X — J (Z (t) : < i < r„) / {Tm < oo) 

Po [Tm < oo) 

(14) = TTTT^ ^ X —7 rl (Tm < oo) . 

Po(T„<oo) Eti^fcexp(^*Z,(r^)) 

Now, clearly on the set {Tm < oo}, there is an index L such that exp {O^^Zl (Tm)) > 

exp(02™')) therefore 

(15) 

'——I (T„ < ») < = ± exp (-«;,n) < 1, 

22k=i '^k exp [eiZk (Tm)) WL .^^ 

where the last inequality follows by Condition 1. Consequently, plugging 

(15) into (14) we obtain that 
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We now are ready to summarize our acceptance / rejection procedure and 
the proof of its validity. 

Algorithm 1.1.1: Simulation of Paths Conditional on Tm < oo 
Step 1: Sample (Z (t) : < t < Tm) according to Pq (•) as indicated via 
equations (11) to (13), and (10). 

Step 2: Given (Z (t) : < t < Tm), simulate a Bernoulli / with probabil- 
ity 

1 

Eti^fcexp {eiZkiTm))' 

(Note that the previous quantity indeed is less than unity due to (15).) 

Step 3: If / = 1, output (Z (t) : < t < Tm) and Stop, otherwise go to 
Step 1. 

As the next result shows, the output of the previous procedure indeed fol- 
lows the distribution of (Z (t) : < t < Tm) given that Tm < oo and Z (0) = 
0. Moreover, the random variable I has exactly probability Pq (Tm < oo). 
So, this procedure actually allows to execute both Steps 4 and 5 in Algo- 
rithm 1.1 simultaneously. We summarize the validity of the procedure in the 
next proposition. 

Proposition 3. The probability of accepting a path in the previous pro- 
cedure is precisely Pq {Tm < oo). In other words P {I = 1) in Step 3 (uncon- 
ditional on the proposed path) of Algorithm 1.1.1 is precisely Pq {Tm < oo). 
Finally, the output of Algorithm 1.1.1 follows the distribution of {'L{t) : < t < Tm) 
conditional on Tm < oo and Z (0) =0. 

Proof. The result follows directly from the theory of acceptance / re- 
jection (see Asmussen and Glynn [2007] page 39-42). In particular, it is well 
known that the bound (16) implies that the probability of accepting a pro- 
posal is indeed Pq (^m < oo). The rest of the statement of the proposal 
follows from our derivation in the construction of Pq (•). □ 

Remark: As mentioned earlier, Assumption A3b) is a strengthening of 
Assumption A3). We can carry out our ideas under Assumption A3) as fol- 
lows. First, instead of (M {t) : t > 0), given a vector a = (ai, 02, a^)"^ with 
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non-negative components that we will explain how to choose momentarily, 
consider the process Ma (•) defined by 

Ma (t) = max (Z (s) + as) . 

s>t 

Note that we can simulate (M (t) : t > 0) jointly with (Z {t) : t > 0) if we are 
able to simulate (Ma (t) : t > 0) jointly with (Z (t) : t > 0). Now, note that 
tfji (•) is strictly convex and that d^i (0) /d9 < so there exists Oj > large 
enough to force the existence of 9* > such that E exp {6*Zi {t) + aiO*t) = 1, 
but at the same time small enough to keep E (Zi (t) + ait) < 0; again, this 
follows by strict convexity of (•) at the origin. So, if Assumption A3b) does 
not hold, but Assumption A3) holds, one can then execute the Algorithm 
1.1 based the process Za (t) := Z (t) + at. 

2.4. Computational Complexity. In this section we provide a complexity 
analysis of our algorithm. We first make some direct observations assuming 
the dimension of the network remains fixed. In particular, we note that 
the expected number of random variables simulated has a finite moment 
generating function in a neighborhood of the origin. 

Theorem 1. Suppose that Al) to A3) are in force. Let r he the coales- 
cence time and N he the numher of random variahles generated to terminate 
the overall procedure to sample Y (oo). Then, there exists 6 > such that 

Eexp {6t + 6N) < oo. 

Proof. This follows directly from classical results about random walks 
(see Gut [2009]). In particular it follows that £'Q(exp (^T^)) < oo. The rest 
of the survey follows from elementary properties of compound geometric 
random variables arising from the acceptance / rejection procedure. □ 

We are more interested, however, in complexity properties as the network 
increases. We shall impose some regularity conditions that allow us to con- 
sider a sequence of systems indexed by the dimension d. We shall grow the 
size of the network in a meaningful way, in particular, we need to make sure 
that the network remains stable as the dimension d increases. Additional 
regularity will also be imposed. 

Assumptions: 
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CI) There exists > independent of d, such that (/ — Q^^ -EX (1) = 

Ad (/-Qj)"'i?Wd(l)-r < -5ol. 

C2) In addition, there exists H, a £ (0, oo) such that 

sup max E[Zf{l)] < a for all i 

d>l l<i<d 

sup max -Eexp {SjWj (1)) < H < oo 
d>i '^<3<d. 

and 3a > independent of d, such that 5j > 9*^ + a for all 1 < j < d 
C3) Finally, C ^ A^ < 1/C uniformly over d>l for some Q G (0,oo). 

As discussed in Section 2.3.1, in Algorithm 1.1, we actually do Step 4 and 
Step 5 simultaneously, therefore, we can rewrite Algorithm 1.1 as follows: 

Algorithm 1.1': Simulate the Coalescence Time 

1. Set r = 0, Z(0) = 0, iV = 0. 

2. Simulate a sample from W — C/r, here U distributed Exp{X) and is 
independent from W. Reset N ^ N + 1, Z{t + U) ^ Z(t) + W-Ur, 

T ^T + U. 

3. If there exists some component i, such that Wi — Uri > —m, return 
to Step 2. 

4. Otherwise, simulate a random walk {C(n)} such that C(0) = and 
C(n) = C(n - 1) + W'(n) - U'{n)r, where W'(n) - U'{n)r are in- 
dependent and identically distributed as W — U'r under the tilted 
measure P' defined in Section 2.3.1. Perform the simulation until 
Nm = inf{n > : Cj(n) > m for some i}. 

5. Reset N ^ N + Nm- Compute p = Y.k=i'^kew{0*kCk{Nm)) and 
sample a Bernoulli / with probability p. If / = 1, Z(t + ^^]^ U'{k)) = 
Z(r) + C{Nm) and r = T + J2k=i U'{k). Return to Step 2. 

6. If / = 0, stop and output r with (Z(t) : < t < r). 

In this algorithm, the total number of random variables required to gen- 
erate is N . Use N{d) instead of N to emphasize the dependence on the 
number of dimensions d. The following result shows that our algorithm has 
polynomial complexity with respect to d: 

Theorem 2. Under assumptions CI) to C3), 
E[N{d)] = 0{d^) asd^oo, 
for some 7 depending on 5o, a, H , a and C,. 
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Denote the number of Bernoulli's generated in Step 5 by Nb and the 
number of random variables generated before executing Step 4 in a single 
iteration by Na- By Wald's identity, we can conclude: 

E[N{d)] = E[Nb]{E[Na] + E[N,n]). 

The following proposition gives an estimate for E[Nm]- 

Proposition 4. Under our assumptions CI) to C3), 

E[Nm] = 0{dlogd) 

and the coefficient in the bound depends only on 6q, a, H and a. 

Proof. First, let us consider the cases in which Wi are uniformly bounded 
from above by some constant B. 

Define ji' = Wi'il}[{9*) / \. Under the change of measure, -E'oEf=i ^'^(n)] 
Sj=i '"^iV'K^i*)/'^ ~ "-A*'' hence Ylii=i Ci{n) — n//' forms a martingale. Nm 
is a stopping time and Ci{Nm) < m + B. By Optional Sampling Theorem, 
we have 

E[N^] = < ^ . 

Using Taylor's expansion, we have 

Un) = MO) + oMio) + ^V'i'lO, 

for some u G [0, 1]. As ^i{0*) = tpi{0) = 1, we have, 

«(o)+^V':'«*)=o. 

As e* > 0, 

(17) ^^(o) + ^V'i'«) = 0. 

Under Assumption CI), ^^(0) = E{Zi{l)) < —Sq. Under Assumption C2), 
we have 

1. i;i{6^) = E[exp{6iZ,{l))] < i/ < oo; 

2. Si > 9* -\- a for all component i and some fixed a > 0. 
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As a result, 

i^';{uet) = E[Z,{lf eMnO*Z,{l))] 

< E[Zi{lfl{Zi{l) < 0)] + E[Ziil)^ exp{9*Z,il))I{Z,iO) > 0)] 

< E[Zi{lf]+E[Z,(lfexp{e*Zi{l))I{Zi{0) > 0)]. 

Note that 

E[Zi{lfeMO:Z,{l))I{Z,{0) > 0)] = S[Z,(l)2exp(-aZ,(l))-exp(,5,Z,(l))(Z,(0) > 0)], 
and for any x > 0, x^exp(— ax) < 4exp(— 2)/q. Therefore, 

^p'-iue*) < E[Z,{lf] + -e-^E[eMS^Zi{m 

a 

a 

Plug this result into equation (17), we get 

(18) ^: > — ra?- 

On the other hand, by a Taylor's expansion again, 

for some u G [0, 1], concluding 

i;-{ue*) = E[Z,il)^expiue*Z,im > E[Z,{lfl{Z,{l) < 0)] 
>i?[r2/(^>l)]=r2exp(-A). 

Therefore, we have 
Since Wi >0 and J2i '^i = 

oo mmj r ■ e 



Now, we can conclude that 

i^[a(n)] (i(m + i3) Arf(m + i3)(fl + ie-^g) 

-C/^VrnJ = < ; < 2 X • 

^ fi do mmj rfe ^ 
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By Assumption C3), we have A < l/Q. Combining Assumption CI) and 
the fact that (/ - Q^)"^ > 0, A > and E[Wi] > 0, we have n > 6o. 
Therefore, we can conchide that E[Nm] = 0{dlogd). Here the term logd 
comes from the fact that m = 0(log(i/minj 0?) as discussed in Section 2.3 
and that 0*'s are positive and uniformly bounded away from as shown by 
(18). 

Now, let's consider the more general cases when the VFj's are not bounded 
from above. For any B > 0, define the truncation Wi = WiI{Wi < B), let 
W* be Wi under exponential tilting by 0* applied to Wi, and define the 
random walk Ci{n + 1) = Ci{n) + Wj(n)*_^_]^ — U*{n)ri. Let = inf{?i : 
Ci{n) > m for some i}. Since C.j(n) < Ci{n), we have Nm < Nm- Since W*^s 
are bounded from above by B, applying the previous result, we have: 

E[N^] < , + . . 

Eti w^E[{Wi - Un) ^Mni^i) - Un)] 

Since —Uri > 0, we have 

E[{Wi-Uri)eMOl{W^-Uri)] 
< E[{W, - U,r) eMG*i(Wi - Un)] - E[Wi eMO*iWi)I{Wi > B)], 

and for all i? > l/a, 

E[WieMO*Wi)\w,>^)\ <Bei^^{-aB)E[e^^{5iWi)] < Bexp{-aB)H. 

Recall that 

V- , ,,,,,,, ,^ Sorfe ^ 



: w,E[{W, - U,r) eMOKW^ - Un)] = ^' > 



^ . . . ....... , + 

Therefore, we can take B = O ( — ^ log ( /"""4^'_!2 h-n ) ) which is independent 
of d and get 

, , ~ , 2Xd(m + B)(a + ^e~^H) 
E[N„,] < E[Nm] < \ /\ \ = 0{d\ogd). 

□ 



Now, we give the proof of the main result in this subsection. 
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Proof to Theorem 2. Recall that 

E[N]=E[B]{E[Na] + E[N^]). 

Since B is the number of trials required to obtain / = 0, E[B] = 1/P{I = 0). 
As discussed in section 2.3.1, P{I = 1) > 1 — J2i=i exp(— therefore 

E[B] < . 

l-Etiexp(-e» 

Similarly, we have E[Na] = 1/P{U > (m + Wi)/ri, V i). For any K > 0, 

> V ^) > PiU > ^^■,W, < K). 

ri miUj r j 

Under Assumption C2), we have 

d 

P{Wi <K)>l-Y^ P{W^ > K)>1- dHexp {-Ka). 

1=1 

As U and W are independent, 

P{U > —— — V f) > exp (m + K)]{1- dexp{-Ka)). 

ri V Cmmjri / 

As discussed in section 2.3.1, we take m = 2 log d/ miuj 0? so that E[B] < 
1/(1 - 1/d). Also, picking K = 2(log(i + logH)/a, we have P{U > (m + 
Wi)/ri, V > (1 - l/(i)exp(-log(i/(Cminirimini6'*)) and hence E[Na] < 
(1 + 1 /d) exp (log d/{C miuj rj miuj 9*)). By the previous proposition we have 
E[Nra] = 0{dlogd). In summary, we have: 

1 1 - 1 

E[N] = E[B]{E[Na]+E[N„^]) < (l+-)(d'^™i ™' +0((i)) = o (d"^ ). 

As discussed in the Proof of Proposition 4, 9* > So/{a + 4exp{—2)H/a) 
and Tj > 1 /(5o are uniformly bounded away from 0, therefore. 



E[N] = 0{d^+^ 

□ 
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3. Algorithm for Reflected Brownian Motion. In this section, we 
revise our algorithm and explain how we can apply it to the case of reflected 
Brownian motion. 

In this setting, we let 

X (t) = vt + AB (t) , 

where v G is a drift vector and A ■ = S € W^^"^ is a positive definite 
covariance matrix and consider the solution to the Skorokhod problem 

dY {t) = dX {t) + ML {t) , Y (0) = yo 
Y {t) > 0, Yj (t) dLj (t) > 0, Lj (0) = 0, dLj (t) > 0. 

We assume that R is an M-matrix, that is, R = I — Q, where Q has non- 
negative components and a spectral radius less than unity, so R~^ has only 
non- negative elements and we also assume the stability condition i?^^v < 0. 
As discussed in the Harrison and Rieman [1981] there is a unique solution 
pair (Y, L) to the Skorokhod problem associated with X, and the process 
Y is called a reflected Brownian Motion (RBM). We wish to sample Y (oo) 
(at least approximately, with a pre-defined controlled error). 

We attempt a construction based on Algorithm 1 in a way that is com- 
pletely parallel to the case in which X (•) is a compound Poisson process. By 
Lemma 1 we can find z such that v < z and R~^z < 0, and we can consider 
the process 

Z (t) = X (t) - zt := fit + AB (t) , 

where ^ = v — z. It follows also from Lemma 1 that if one considers the 
Skorokhod problem (4), which has orthogonal reflection, then R~^Y (t) < 
R~^Y+ (t). 

Now, in order to apply the strategy used in the previous sections to RBM, 
we need to address two problems. First, the input process Z requires a con- 
tinuous path description while the computer can only encode and generate 
discrete objects. Second, the dominating process is a reflected Brownian mo- 
tion with orthogonal reflection, therefore the hitting time r to the origin is 
infinity almost surely (see Varadhan and Williams [1985]), which means that 
Algorithm 1 will not terminate in finite time in this case. To solve the first 
problem, we take advantage of a wavelet representation of Brownian mo- 
tion and use it to simulate a piecewise linear approximation with uniformly 
small (deterministic) error. To solve the second problem, we choose r as the 
first passage time to a small ball around the origin and we find bounds on 
E[t] and we also control the error in the re-construction of the Skorokhod 
problem forward in time. In sum, we concede to an algorithm that is not 
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exact but one that could give any user-defined e precision. Nevertheless, at 
the end of Section 3.1 we will show that we can actually use this e-biased 
algorithm to estimate without any bias the steady-state expectation of any 
continuous function of RBM by introducing an extra randomization step. 

Section 3 is organized as follows. In Section 3.1, we will describe the 
main strategy of our algorithm. In Section 3.2, we will develop an algorithm 
to simulate a piecewise linear approximation of Brownian motion using a 
wavelet representation. In Section 3.3, we will discuss the details in simu- 
lating the stationary dominating process combing a wavelet representation 
with the techniques we have already used for the compound Poisson cases. 
In the end, in Section 3.4, we will give an estimate of the computational 
complexity of our algorithm. 

3.1. The Structure of the Main Simulation Procedure. The main strat- 
egy of the algorithm is almost the same as Algorithm 1 except for two 
modifications due the two issues discussed above: first, instead of simulating 
the input process Z exactly, we simulate a piecewise linear approximation 
such that \Zi{t) — Zi{t)\ < e for all components i and t > 0; second, instead 
of sampling the coalescence time r such that M(t) = Z(t), we simulate an 
approximation coalescence time, r^, such that M(re) < Z(re) -|- e. 

With these notations, we now give the structure of our algorithm. The 
details will be given later: 

Algorithm 2: Sampling with Controlled Error of Y (oo) 
Step 1: Let > be any time for which M (t^) < Z (re)-|-e and simulate, 
jointly with r^, Z^^^ (t) = -Z" (r^ - t) for < t < Te. 

Step 2: Define Xt;^ (t) = -Z^ (r, - t)-zx (r^ - t) and compute Yl^^ (r^, 0) 
which is obtained by evolving the solution Yl^^ ("i 0) to the Skorokhod prob- 
lem 

dYt^^ {t, 0) = dXr^^ (t) + RdL^r {t, 0) 
Y%^ (t, 0) > 0, Yl.^j (t, 0) dL^r.j (t, 0) > 0, L^r.j (0, 0) = 0, dL_^^,j {t, 0) > 0, 

for Te units of time. 

Step 3: Output Y1^^(t,,0). 

First, we show that there exists a stationary version {Y*(i) : t < 0} that 
is coupled with the dominating stationary process {Y"'"(t) : t < 0}. 

Lemma 2. There a stationary version {Y*(t) : t < 0} of Y such that 
R-^Y*{t) < R~^Y+{t) for all t < 0. 
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Proof. Using similar notations as in Section 2.1. Let Yxi-jy) be the 
solution to the Skorokhod problem: 

dYT {t, y) = dZ^ (t) + RdLr {t, y) , Yt (0, y) = y, 
Yt (t,y) > 0, Ytj (i,y) dLrj (t,y) > 0, Ltj (0,y) = 0, dLrj (t,y) > 0. 

By Lemma 1, we have i?-^Yr(t, Y+(r)) < R-^Y+{T + t) for ah < t < 
|T|. Moreover, for any Ti > T2 > t > 0, we have 

ii-^Y_T2(T2 -t,Y+(-r2)) > i?-iY_Ti(T2 -t,Y„Ti(ri -r2,Y+(-Ti))) 

= i?-iY_Ti(Ti-t,Y+(-ri)) 

Define Yrit) = Yrit - T, Y+(-T)) for ah T < t < 0. From the above 
inequality, we know that for fixed t, R~^YT{t) decreases as T — t- —00. On 
the other hand, by the property of Skorokhod problem, R~^YT{t) > 0, 
therefore, lim7^^_oo i?~^Yj'(t) is well defined in M"^. As R invertible, we can 
define Y*(t) = liuiT^-ooYTit). Using the same argument as for DCFTP in 
Kendall [2004], we can conclude that Y*{t) is a stationary version of Y and 
R-'^Y*{t) < R-^Y+{t). □ 

The following proposition shows that the error of the above algorithm has 
a small and deterministic bound. 

Proposition 5. Suppose X £ R''. Let r = maxij R~j / mmij{R^j : 
R^j > 0}. Then there exists a stationary version Y*ofY such that in each 
component i 

\Y*{0)-Yl,{n,0)\<{l + dr)e. 

Proof. Consider three processes on [— r^, 0]. The first is the coupled sta- 
tionary process Y*(-) as constructed in Lemma 2, which is the solution to 
the Skorokhod problem with initial value Y*(— r^) at time — and input 
process X(-) = X(Te) — X(— •) on [— r^, 0]; the second is a process Y(-), which 
is the solution to the Skorokhod problem with initial value at time — and 
input process X(-); the third is the process Yl^^(t, 0) as we described in the 
algorithm, which is the solution to the Skorokhod problem with initial value 
at time — and input process X^^^ (t) as defined in Step 2 of Algorithm 
2. 

By definition, we know that for each component i, \Y-'^{—Te)\ < e. Since 
i?"^Y(T,) < R-'^Y+iTe), the coupled process Y*{-Te) < dre. Note that 
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Y*(-) has the same input data as Y(-) except for their initial values. Ac- 
cording to the comparison theorem of Ramasubramanian [2000], the dif- 
ference between these two processes is uniformly bounded by the differ- 
ence of their initial values in each component. Therefore, we can conclude 
|y,*(0)-y,(0)| <dre. ^ 

On the other hand, Y(-) and Yl^^(-,0) have common initial value and 
input processes whose difference is uniformly bounded by e. Since the Sko- 
rokhod map is a contraction under the topology of uniform convergence on 
compact sets, we can conclude that \Yi{0) — Y^^ ^(r£,0)| < e. 

Simply apply the triangle inequality, we obtain that 

\Y*{0)-Yli{n,0)\<{l + dr)e. 

□ 



We conclude this subsection explaining how to remove the e-bias induced 
by Algorithm 2. Let T be any positive random variable with positive density 
{f{t) : t > 0} independent of Y*(0). Let g( : R"^ — >• R be any positive function 
such that \g{x) — g{y)\ < K maxi=i \xi — yi\ for some i^' > 0. As illustrated 
in Beskos et al. [2012], 





r9(y*m 






E 


/ dt 


= E 






Jo 




Jo 



m 
fit) 



dt 



l(g(Y*(0)) >r) 
f{T) 

Since |5^j*(0) — Y^^-{t^,0)\ < (1 + dr)e, we can sample T first, and then 
select e > small enough, output l(g'(Y;^^ (r^, 0)) > T)/f{T) as an unbiased 
estimator of E[g{Y*{0))] without the need for computing Y*(0) exactly. It 
is important to have (Y^^{t^, 0)) coupled for e — )■ and this can be achieved 
thanks to the wavelet construction that we will discuss next. 



3.2. Wavelet Representation of Brownian Motion. In this part, we give 
an algorithm to generate piece-wise linear approximations to a Brownian 
motion path by path, with uniform precision on any finite time interval. 
The main idea is the use of a wavelet representation for Brownian motion. 

By the Cholesky decomposition, any multidimensional Brownian motion 
can be decomposed as linear combination of independent one-dimensional 
Brownian motions. Now, suppose we want to give a piece- wise linear approx- 
imation to a d-dimensional Brownian motion Z with uniform precision e on 
[0, 1]. Suppose we can write Z = AB, where A is the Cholesky decomposition 
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of the covariance matrix and the Sj's are independent standard Brownian 
motions. If we are able to give a piece-wise hnear approximation Bi to each 
Bi on [0,1] with precision e/{d ■ a), where a = maxjj \Aij\, then AB is a 
piece-wise Hnear approximation to Z with uniform error e. Therefore, in the 
rest of this part, we only need to work with a standard one-dimensional 
Brownian motion. 

Now let's introduce the precise statement of a wavelet representation of 
Brownian motion, see Steele [2001] page 34-39. First we need to define step 
function H{-) on [0, 1] by 



Hit) 



1 for < t < i 
-1 for i < t < 1 
otherwise. 



Then define a family of functions 

Hk{t) = 2~^/^H{2H - I) 

for n = 2^ + 1 where j > and < / < 2^'. Set Ho{t) = 1. The following 
wavelet representation theorem can be seen in Steele [2001]: 

Theorem 3. If {W^ : < k < 00} is a sequence of independent stan- 
dard normal random variables, then the series defined by 

Bt = Y.[^l Hk{s)ds 

k=0 ^ 







converges uniformly on [0, 1] with probability one. Moreover, the process {Bt} 
defined by the limit is a standard Brownian motion on [0, 1]. 

Choose % = 4 • ^/log k and note that -POVF'^'J > %) = 0{k~'^), so 
Efclo^d^'^l > Vn) < 00. Therefore, P{\W''\ > ?7fc,i.o.) = 0. The sim- 
ulation strategy will be to sample the W'^^s jointly with the times, k, at 
which ]VK^] > 

Note that if we take j = [log2 k], as shown in Steele [2001], 



(w' / H,{s) ds) < (2-^/^ . 



max \W'^ 

2:'<fc<2J + l-l 



Since '^j=o 2 -'^^Vi + 1 < 00, for any e > there exists Kq > 0, such that 
Ej = riogXol 2~^/2vJ+T < e. As a result, define K = max{k : \W''\ > 
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Vk} V Ko < oo, then Y2'k=K+i 

we can simulate 

{iW')tvI<} jointly, 

(19) B^{t) = y^W^ / Hk{s)ds 

fc=o 

will be a piecewise linear approximation to a standard Brownian motion 
within precision e in C[0, 1]. 



Define 1^ as a Bernoulli with probability p{k) = P{\W^\ > rji^), then 
= W^Ik + W^{1 - Ik), where W'' follows the distribution of W'' given 
IW'^I > rjk and W'' follows the distribution of given \W^\ < rj^. Define 
Go = and = inf{A; > Gi : Ik = 1}, then K = max{G/ : Gi < oo, / > 
0} < oo. Once Kq = Gq < Gi < ... < < = oo are simulated, we can 
then generate {VF'"'} up to K = by: 




W'' If n = Gj for some j < I — I; 
Otherwise. 



Then we can evaluate B'^{t) according to (19), and \B{t) — -6*^(^)1 < e for all 
tG [0,1]. 



Now we describe the procedure to simulate Kq = Gq < Gi < ... < 
G;_i < Gi = oo. Since the W^^s are independent, we only need to describe 
the procedure to simulate Gi which is basically the same as that of Gk given 

Gk-l- 

Note that P(Gi = oo) = \{t=Ko+ii'^ " P'^^)) ^ [0' 1] p{k) < fc"^, 
therefore, for any > 0, we have: 

N N oo 

n {i-Pik))>p{G^ = ^)> n a-pimi- E ^~') 

k=Ko+l k=Ko+l k=N+l 

N 

>(l-iV-3) -Q (l_p(A;)). 

fc=-fs:o+i 

To simulate a Bernoulli with P{Gi = oo), we simply simulate 

/(^7 < P(Gi = oo)), 

where U is uniformly distributed on [0, 1] and the condition U < P{Gi = oo) 
can be obtained by sending N large enough without computing the infinite 
product \{t=Ko+ii'^-P^k))- 
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If U > 0^^^-0+1(1 we need to simulate the distribution with 

probabihty mass function: 

P(C. - i\G. < ool - ^(^')nCk+i(l-P(^)) < .-4 . ^ 
P(Gi-j|Gi<oo)- ^_^oo^^(^_^(^)) <J A, 

where A = — Y[T=Ko+i(^ is a finite constant. As a consequence, 

we can use acceptance/rejection to simulate Gi. In detail, we first propose 
G > Kq + 1 with probability mass function proportional to and then 
compare U, a uniform on [0,1] independent of G, with the likelihood ratio: 

p(g)nt-xo+i(i-pW) 

LK[G) = . 

We repeat proposing G and sampling U until U < LR{G) at which point 
we output Gi = G. 

If [/ < (1 - N-^)Y{^^^^^^{1 - p{k)), we conclude d = 00 and stop. 



Also we note that B''{t) has the following nice property: 



Proposition 6. 

B\\) = B{1). 

Proof. The equality follows immediately from the fact that Hn{s)ds = 
for any n > 1 and m > 1. □ 

As a consequence of this property, for any compact time interval [0,T], 
(without loss of generality, assume T is an integer), in order to give an 
approximation for B{t) on [0,T] with guaranteed e precision uniformly in 
[0,T], we only need the run the above algorithm T times to get T i.i.d. 
sample paths {B'^''^'^\t) : t £ [0, 1]} for i = 1,2, ...,T, and define recursively: 

B^it) = ^B^^^^\l) + Bl^it-[t\). 
1=1 

3.3. Joint Simulation of and . Our goal now is to develop an algo- 
rithm for simulating and (Z'^(t) : < t < r^) jointly. In other words, we 
simulate Z'^(t) forwards and try to find a random time such that for any 
time s after r^, < Ziir^ + e. 
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Now suppose {W^{i) : n,A; G N, i = 1,2,..., d} are i.i.d. standard normal 
random variables. By the wavelet representation we can write, for any t = 
n + s, s £[0, 1]: 




Let us put (20) in matrix form: 

oo s 

(21) Z(t) = Z(n) + ylW° -s + ^VwJ;- / Hk{u)du 

For all n > 0, let B„(s) = A Y.t=i !o Hk{u)du, then {B„(s) : n > 0} are 
i.i.d. linear combinations of Brownian bridges. Note that {W'^{i) : n > 0} 
is independent of {W!^^{i) : A; > 1}. We can split the simulation into two 
independent parts: 

1. Simulate the discrete-time random walk {Z(n)} using {W^{j)}. That 
is Zm = and Zi{n + 1) = Zi{n) + AijW^^^{j). 

2. For each n, simulate B„(s) to do bridging between Z(n) and Z(n+ 1). 

Intuitively, time t has a good chance to be an approximate coalescence 
time if: 1). (Z(n) : n > \t\) never exceeds level Z(Te) + e; 2). B„(s) for n> t 
do not reaching a substantially high value. We shall translate this intuition 
into rigorous mathematics and develop an algorithm to simulate jointly 
with (Z'(t) : < t < Te). 

As we have assumed in Section 2, fii < —5o for some 5o > 0. For conve- 
nience in some calculations that follow, we define C = so in particular 
fii < —5(). Define S(n) = Z(n) + nd^l and note that (S(n) : n > 0) is a ran- 
dom walk with strictly negative drift. Choose m > satisfying Condition 1 
for the random walk (S(n) : n > 0) as discussed in Section 2.3. Define 

(22) Tm = inf{n : 5j(n) > m for some component i}. 

Recall that (B.„(s) : n > O) are the i.i.d. linear combinations of Brownian 
bridges. Define 

(23) M = max{n > : max (B„^j(s) — nQ > for some component i}. 

0<s<l ' 

We will show in Section 3.3.1 that M < oo. 
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If Tm = oo then for all t > M and every component i, 
(24) Zi{t) <m + Bn,i{t-n)-tC<m, 

where n = \ t\. 

Now we can give the main conceptual trust behind our algorithm: 
Algorithm 2.1: Simulating and Z*^ Jointly. 

The output of this algorithm is (Z^ : < t < r^), and the approximation 
coalescence time r^. 

1. Set = 0. 

2. Simulate Z*^ forwards until some integer time ni > r,. such that there 
exists no < ni satisfying Zl{s) < Z^^uq) + e for all no < s < ni and 
Zi{ni) < Zi{no) — m for all i = 1, d. Set ^ uq. 

3. Simulate a Bernoulli Ii with parameter p = Po(^m < oo), Tm. as 
defined in (22). 

4. If /i = 1, simulate a new conditional path (C(n) : < n < T^) 
following the conditional distribution of {S{n) : < n < Tm) given 
that < oo and S(0) = 0. Let Z'(ni + n) = Z^(ni) + S(n) - nC for 
< n < Tm. For all t G [ni, ni -|-Tm], simulate Z'^(t) according to (20). 
Set Te ni + Tm and return to Step 2. 

5. If /i = 0, simulate a new sequence of bridges (B„(s) : < n < M), 
where M is defined in (23). Simulate (S(n) : < n < M) following the 
conditional distribution given that Tm = oo and S(0) = 0. Compute 
Z%t) for all t G [ni,ni + M] according to (20). 

6. If there exists < I < M and some component i such that Zi{ni+l) > 
Zi{no) + e, set ^ ni + I and return to Step 2. 

7. Otherwise, we have Zi{t) < Zj(to) + e for all t G [ni,ni + M]. On 
the other hand, according to (24), Zi{t) < Zi{nQ) for all t > ni + M. 
Therefore, we can output = no as the approximate coalescence time 
jointly with (Z'(t) : < t < rj. 

Since S(l) is a multidimensional Gaussian random vector. Assumption 
A3) and Condition 1 are satisfied. Therefore, we can use the same tech- 
nique to implement Step 3 and Step 4 as in Algorithm 1.1. The rest of this 
subsection will be focused on Steps 5 and 6. 

3.3.1. Simulating M and (B^(i) : < n < M). Recall that B„(t) = 
^Z^fe^i^n ■ Jq Hk{u)du, where W^{i) are i.i.d. standard normal for each 
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component i. Note that 

E E ^(i^n «i > 4\/i^ + 47!^) < E E 7^ < ^- 

n k n k ^ ' 

By the Borel-Cantehi Lemma, we can conclude that for each i G {1, d} 
there exists an integer random variable M* > such that for all nk > M*, 
|l^n(OI ^ 4T/log n + 4\/log /c. Clearly, i/log t = o(t) as t — )• oo, so we can 
select a mo large enough such that for any n > rriQ 

oo 

nC - a(i(4v^logn - E 2"-' V^) > 0. 
Then, for any n > max^^^^ M* V mo, 

oo „f 

B„(t) = ^EWn- / ^fc(n)dn 

oo 

< ad(4yi^ + E 2"^/'\/i) < nC, 
i=i 

where, j = [log2 A;] . Therefore, we can choose M = maxj M* V mo- Based 
on this, we can now give the main procedure to simulate and (Z(t) : < 
t < Te) jointly: 

Algorithm 2.1.1: Simulating of M and (B^(i) : < n < M) jointly 

1. For each component i, simulate Mi and : n > 0,k > l,nk < 
M). Compute M = max, M* V mo 

2. For each < n < M and each component i simulate = max{A; > 
Kq : W!^{i) > 4-v/log A;}, where JCo is defined as in Section 3.2, and 
(VF^(i) : 1 <k < K^), where W^(i) for k < Mi/n are already given in 
Step 1 while W^ii) for k > Mi/n are independent standard normals 
conditional on {|VF'^(i)| < 4(^/logn + Vlog A)}. 

3. For any < n < M, compute and output 

^n,(i) = E^^i(E^"« / Hk{u)du). 

i=i k=i 

To execute Steps 5 and 6 in Algorithm 2.1, we first use Algorithm 2.1.1 
to simulate M and (B^(t) : < n < M). As discussed previously, the 
Brownian bridges B^(t) are independent of the random walk S(n). For M 
given, if we can simulate a path (S(n) : < n < M) conditional on = oo 
and S(0) = 0, Z^(t) can be computed according to (20). 
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3.3.2. Keeping Track of the Conditioning Events. Algorithm 2.1 requires 
keeping track of various conditioning events in two cases. First, when Ii = 0, 
we need to simulate S(n) conditional on that Tm = oo as just mentioned. 
Second, to simulate the Brownian bridges, we might need to sample W^{i) 
from a truncated standard normal instead of a standard normal, as described 
in Step 2 of Algorithm 2.1.1. 

In order to deal with the first case, we use the acceptance-rejection ap- 
proach given in Blanchet and Sigman [2011] by which we can simulate {S(n)} 
jointly with two sequence of random times {Ti : / > 1} and {A; : / > 1} 
defined as following: 

1. Set Ai = min{n : Si{n) < —2m for every i}. 

2. Define Ti = min{n > A/ : Si{n) > Si{Ai) + m for some i}. 

3. Put A/+1 = min{n > TiI{Ti < cx3)vA/ : Si{n) < Si{Ai)—2m for every i}. 



-2m 



S{A2) + m 



■ 2m ' 



-5- 
-6 



Ai 



Fi 



Ao 



10 11 



12 



As illustrated by this picture, at t = A/, /i = if F/ = oo and /i = 1 
if F; < oo. Moreover, (S(n + A;) - S(A/) : n > 0) for F; = oo follows 
the conditional distribution of (S(n) : n > 0) given that S(0) = and the 
random walk never exceeds level m, see Lemma 3 in Blanchet and Sigman 
[2011]. The approach of Blanchet and Sigman [2011], which works in one 
dimension, must be modified using the change-of-measure as described in 
Section 2.3.1. 

For the second case, we just need to keep track of the bounds for each 
Wn{i) throughout the algorithm. Now, we can write down the precise ver- 
sion of Algorithm 2.1: 
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Algorithm 2.1': Simulating and (Z^(t) : < t < r^). 

The output of this algorithm is (Z*^ : < t < r,.), and the approximation 
coalescence time r^. 

1. Set B^{i) = oo for all n > 1, /c > 1 and 1 < i < (i. Set L = and 
r, = 0. 

2. Simulate S(n) until A/, where I = min{j : Tj = oo, Aj > t^}. Compute 
Z^(n) = S(n) - nC 

3. Do bridging for t G [t^, A/] given < B^^{i) and compute Z'^(i) 
according to (20). If there exist some t > r;_i such that for all t < 
s < A;, Zf{t) > Zf{s) - 2e and > (A;) + m - 2e, set r, ^ t 
and go to Step 4. Otherwise, set A/ and return to Step 2. 

4. Simulate M jointly with given < Set ^ 
4Vlog n + CVlog k for ah n • A; > M*. 

5. Do bridging for t G [A;, A; + M] given < B^{i) and compute 
Z'^(f). If there exist some t and i such that > Z^{t^) + e, set 
Te ^ t and return to Step 2. 

6. Otherwise, stop and output as the approximation coalescence time 
along with (Z^(t) : < t < r,). 

3.4. Computational Complexity. In this part, we will discuss the com- 
plexity of our algorithm when d and the other parameters /i and A are fixed 
but send the precision parameter e to 0. Denote the total number of random 
variables needed by A^(e) when the precision parameter for the algorithm is 
e. 

Recah that Z(t) = -^t + AB{t). We set 5, a G (0,oo) such that: 
Dl) > 5 > 0; 
D2) maxj Aij < a. 

The following result shows that our algorithm is polynomial in e: 
Theorem 4. Under Assumptions Dl) and D2), 

E[N{e)] = 0{e--^-^log{-^)), as e ^ 0, 
where ac is a computable constant depending only on A. 

The random variables we need to simulate in the algorithm can be divided 
into two parts: 
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1) The random variables used to construct the discrete random walk Z(n) 
for n < T^; 

2) The conditional normals used to bridging between Z(n — 1) and Z(n). 

Since > ij) and < B) are negative correlated, it follows that 

P{\W\ > r]\ \W\ <B)< P{\W\ > ??). 

Therefore, the expected number of conditional Gaussian random variables 
used for Brownian bridges between Z(n — 1) and Z(n) is smaller than the 
expected number that we would obtain if we use standard Gaussian random 
variables instead. 

Let K = max{A; : \Wk\ > %} V Kq as defined in Section 3.2. As discussed 
above, the expected number of truncated Gaussian random variables needed 
in part 2) is bounded by 

Therefore, we conclude that 

E[N{e)] < idE[K] + l){E[T,] + l). 
To prove Theorem 2, we first need to study E[K] and E[Tf]. 

Proposition 7. 

E[K] = 0{e-^ log (^^^) 

Proof. Recall that r]k = 4A/log k, let pk = P{\W''\ > r/^), then pk = 
0(/c~^). Therefore 

oo oo oo 

E[K] = ^P{K >n)<Ko+ ^ 

n=l n=Ko+l k=n 

oo oo 

= Ko+ Yl k-pk<Ko + 0{Y,k-''). 

k=Ko+l k=l 

The second term of the left-hand side is finite and independent of e and Kq . 
On the other side, 

i=log2 Ko ^ ^ 

Therefore, we can choose Kq = 0{e~^ log (^)) such that ^j^y^g^ Kq 2~-'^^Vj + 1< 
e. 
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In order to get the approximation within error at most e for the d di- 
mensional process, according to the Cholesky Decomposition as discussed in 
Section 3.2, we should replace e by J^. Therefore, 

m] = 0{(4-Y"^oJ^)) = 0(e-^ log (^)). 



da J \ e / \ e , 

□ 



The rest we need to do is to estimate Tg.Let be the time before the 
algorithm executes Step 4 in a single iteration. Using the same notation as 
in Algorithm 2.1 and a similar argument as in Section 2.4, we have 

. , ^ E[Ta] + E[T^\T^ < oo] + E[M] 
^^'^ P{Tm < oo)p 

Where 

p = P(maxZ|(t) <m + e,VO<t <M |Z(0) =0;S{n) <m). 

i 

It follows easily that 

p > P(maxmaxZi(t) < m|Z(0) = 0). 

i t>0 

Since S(l) is a multidimensional Gaussian random vector, Assumptions 
C2) is satisfied, along with Assumption Dl) and D2). Applying Proposi- 
tion 4, we can get upper bounds for E[Tm\Tm < oo], l/P{Tm < oo) and 
1/P(maxj maxt < m|Z(0) = 0), which depend only on d,a and 6 and 
thus are independent of e. Besides, the bond on E[M] can be estimated by 
the same method as in Proposition 7 in terms of C = <^/2, hence is also 
independent of e. Therefore, we only need to estimate -E[Ta]. 

Proposition 8. E{Ta\ = 0{e~"'^) as e — > 0. Here ac only depends on 
the matrix A. Moreover, in the special cases where Aij > 0, ac = d. 

Proof. RecaU that Z{t) = -fit + AB{t). We divided the path of Z{t) 
into segments with length ^in^iA- 

r/^/, 2(m + e) , 2(?Ti + e), , 
{{Z{k ■ ^ 7 + s) : < s < -) :k>0}. 
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Let 

Nb = mm{k : AB{k ■ ^(!!L±£l + s) < e for all < s < ^i!!^±£l}. 



By independence and stationarity of the increments of Brownian motion, 
is a geometric random variable (hence is finite) with parameter: 

p = p(ABik.^-^^ + s)<eior ah < . < 

\ 

On the other hand, since /ij < —6, we have 

1. Zm ■ '-^ + s) < Zm ■ + for all < s < 

2. Z,{{N, + 1) • '-^) < Zm . '-^) - m. 

Therefore, Algorithm 2 should execute Step 4 after at most ^(™+'^) ^jy^j^i'^ 
units of time in a single iteration: 

dp 
From this inequality, it is then sufficient to show that p = 0(e"'^). 

Note that the set C = {y £ : Ay < e} forms a cone with vertex A e 
in since A is of full rank under our assumption. Define tq = inf{t > : 
B{t) ^ C} given B(0) = 0, then 

2(m + e), 
p = P{tc > ^ ^ - 

If d = 2, it is proved by Burkholder [1977] that ac = f where 9 G [0, vr) is 
the angle formed by the column vectors of A~^. Therefore, we can compute 
that 

AuA21 + ^12^22 



cos 6 



^{Ai, + Ai,){Ai,+Ai,y 



which only depends on A. 

On the other hand, if d > 3, applying the results on exit times for Brow- 
nian motions given by Corollary 1.3 in DeBlassie [1987], 

P{rc>^-^^)^u.\\A-h\r 



as e — >• 0, here || • || represent the Euclidian norm and u is some constant 
independent of e. The rate ac is determined by the principal eigenvalue of 
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the Laplace-Beltrami operator on (S'^^^ n C), where S''^^ is a unit sphere 
centered at the vertex of C, namely A~^e. The principal eigenvalue only 
depends on the geometric features of C and it is independent of e, hence so 
is ac- Since A is given, we have 

P^rc > = 0(e»^) as e ^ 0. 



For d > 3, computing of ac is not straightforward in general. However, 
when Aij > 0, we can estimate ac from first principles. Indeed, if Aij > 0, 
we have that 

C = {y G : Ay < e} C {y G M'^ : < ^j. 

ad 

As the components of B(t) are independent, 

p>P\ max B(t) <—,] , 

\o<f<^(Y'' / 

where B{-) is a standard Brownian motion on real line. 
Applying the Reflection Principle, we have 

p( max B{t) <±]= ^ =exp ( ) = 0(e). 

As a result, p = 0{e'^) when the correlations are all non-negative. □ 
Given these propositions, we can now prove the main result in this part. 



Proof of Theorem 4. As we have discussed, 

E[N{e)] < {dE[K] + l){E[T,] + l). 
Also by Proposition 7, 

As discussed above, 

E[Ta] + E[T„,\T^ < oo] + E[M] 



E[re] < 



P{Tm < oo)P(maxj maxt>o .^j(i) < m|Z(0) = 0) 
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According to Proposition 8, E\Ta\ = 0{e~°''^) and ac is a constant when A 
is fixed. As we have discussed, E[rm\Tm < oo], P{Tm < oo), P(maxj maxt Zi{t) < 
m|Z(0) = 0) and E[M] are independent of e. Therefore, 



E[T,] = 0{e-'^o- 

In sum, we have 



E[Nie)]=0{e-''<^-Hog Q). 



□ 
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